library(MASS)
library(texreg)
library(pscl)
library(ggpubr)
library(ggeffects)
library(tidyverse)

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

# Load prepared dataset
dat <- read_rds("Data/df_analysis.rds") 

# Regular models ----

#* Geographic references ----
controls <- c("direct", "age2017", "sex", "partylistleader_0", "academic", 
              "renewed_candidacy", "decentralization", "dual_candidacy", 
              "population_density", "household_income", "universityentrancepercent", 
              "youngdemographypercent", "chance_to_win_district_10", "safelist")
model_geoFB <- glm.nb(as.formula(paste0("Posts2constituency ~ ",
                                 paste0(c(controls, "log(fb_freq+1)"), collapse = "+"))),
                      data = dat)
model_geoTW <- glm.nb(as.formula(paste0("Tweets2constituency ~ ",
                                 paste0(c(controls, "log(tw_freq+1)"), collapse = "+"))),
                      data = dat)

#* Regionalized wording ----
model_regFB <- glm.nb(as.formula(paste0("fb_regionalization ~ ", 
                                 paste0(c(controls, "log(fb_freq+1)"), collapse = "+"))),
                      data = dat)
model_regTW <- glm.nb(as.formula(paste0("tw_regionalization ~ ", 
                                 paste0(c(controls, "log(tw_freq+1)"), collapse = "+"))),
                      data = dat)

# Print results ----

htmlreg(list(model_geoFB, model_geoTW, model_regFB, model_regTW),
        custom.header = list("Geographic references" = 1:2,
                             "Regionalized wording" = 3:4),
        custom.model.names = c("Facebook", "Twitter",
                               "Facebook", "Twitter"),
        stars = c(0.001, 0.01, 0.05),
        digits = 2, 
        use.packages = FALSE, 
        custom.note = "Note: Negative binomial regression models with unstandardized coefficients and standard errors. %stars",
        scalebox = 0.9,
        float.pos = "H",
        caption = "Effects of the mandate mode on regionalized wording and geographic references in MPs' social media posts.",
        caption.above = TRUE,
        single.row = TRUE,
        include.adjrs = TRUE, 
        include.nobs = TRUE,
        include.bic = FALSE,
        include.dispersion = FALSE,
        include.groups = FALSE,
        custom.coef.map = 
          list("direct" = "Direct mandate",
               "chance_to_win_district_10" = "Electoral marginality (10% threshold)",
               "safelist" = "Party list electoral viability",
               "decentralization" = "Party emphasis on local politics",
               "CDUCSU" = "CDU/CSU",
               "direct:CDUCSU" = "Direct mandate X CDU/CSU",
               "dual_candidacy" = "Dual candidacy",
               "academic" = "Highly educated",
               "age2017" = "Age (in 2017)",
               "sex" = "Gender (male = 1)",
               "renewed_candidacy" = "Incumbent",
               "partylistleader_0" = "Spitzenkandidat_in (Party list leader)",
               "population_density" = "Population density",
               "household_income" = "Income of private households",
               "universityentrancepercent" = "University entrance qualification",
               "youngdemographypercent" = "Share of young population",
               "log(fb_freq_0 + 1)" = "Facebook frequency (logged)",
               "log(fb_pagelikes_0 + 1)" = "Facebook page likes (logged)",
               "log(fb_comments_0 + 1)" = "Facebook comments (logged)",
               "log(tw_freq_0 + 1)" = "Twitter frequency (logged)",
               "log(tw_followers_0 + 1)" = "Twitter followers (logged)",
               "log(tw_ment_freq_0 + 1)" = "Twitter @-mentions (logged)",
               "log(fb_freq + 1)" = "Facebook frequency (logged)",
               "log(fb_pagelikes + 1)" = "Facebook page likes (logged)",
               "log(fb_comments + 1)" = "Facebook comments (logged)",
               "log(tw_freq + 1)" = "Twitter frequency (logged)",
               "log(tw_followers + 1)" = "Twitter followers (logged)",
               "log(tw_ment_freq + 1)" = "Twitter @-mentions (logged)",
               "(Intercept)" = "Intercept"),
        file = "main_regressionmodels.doc",
        doctype = "doc")


# Prediction plots ----

#* Regular models ----
dat_regFB <- ggeffects::ggpredict(model_regFB, terms = c("direct"), ci.lvl = 0.95) 
dat_regFB$x <- dplyr::recode(dat_regFB$x, "0" = "List MPs", "1" = "Direct MPs")

p.int_regFB <- ggplot(dat_regFB, aes(x, predicted)) +
  geom_point(position = position_dodge(.3), size = 2) +
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high),
    position = position_dodge(.3), width = 0.25) + 
  theme_minimal() +
  ylab('Regionalized wording') + 
  xlab(element_blank()) +
  theme(legend.title = element_blank(), 
        plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  scale_colour_manual(values = c("darkgrey", "black"))


dat_regTW <- ggeffects::ggpredict(model_regTW, terms = c("direct"), ci.lvl = 0.95) 
dat_regTW$x <- dplyr::recode(dat_regTW$x, "0" = "List MPs", "1" = "Direct MPs")

p.int_regTW <- ggplot(dat_regTW, aes(x, predicted)) +
  geom_point(position = position_dodge(.3), size = 2) +
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high),
    position = position_dodge(.3), width = 0.25
  ) + 
  theme_minimal() +
  ylab(element_blank()) + 
  xlab(element_blank()) +
  theme(legend.title = element_blank(), 
        plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  scale_colour_manual(values = c("darkgrey", "black"))


dat_geoFB <- ggeffects::ggpredict(model_geoFB, terms = c("direct"), ci.lvl = 0.95) 
dat_geoFB$x <- dplyr::recode(dat_geoFB$x, "0" = "List MPs", "1" = "Direct MPs")

p.int_geoFB <- ggplot(dat_geoFB, aes(x, predicted)) +
  geom_point(position = position_dodge(.3), size = 2) +
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high),
    position = position_dodge(.3), width = 0.25) + 
  theme_minimal() +
  ggtitle('Facebook') +
  ylab('Geographic references') + 
  xlab(element_blank()) +
  theme(legend.title = element_blank()) + scale_colour_manual(values = c("darkgrey", "black"))


dat_geoTW <- ggeffects::ggpredict(model_geoTW, terms = c("direct"), ci.lvl = 0.95) 
dat_geoTW$x <- dplyr::recode(dat_geoTW$x, "0" = "List MPs", "1" = "Direct MPs")

p.int_geoTW <- ggplot(dat_geoTW, aes(x, predicted)) +
  geom_point(position = position_dodge(.3), size = 2) +
  geom_errorbar(
    aes(ymin = conf.low, ymax = conf.high),
    position = position_dodge(.3), width = 0.25
  ) + 
  theme_minimal() +
  ggtitle('Twitter') + 
  ylab(element_blank()) + 
  xlab(element_blank()) +
  theme(legend.title = element_blank()) + 
  scale_colour_manual(values = c("darkgrey", "black"))

#Plot
theme_set(theme_pubr())
plot.interactions <- ggarrange(p.int_geoFB, 
                               p.int_geoTW,
                               p.int_regFB,
                               p.int_regTW, 
                               common.legend = T,
                               legend = "right",
                               ncol = 2, nrow = 2)
annotate_figure(plot.interactions,
                left = text_grob("Predicted values", rot = 90))
ggsave(file = "effect_plots.pdf", height = 6.5, width = 7)


# Reported predicted effect sizes in Manuscript - regular plots ----

#Facebook Geographic references
dat_geoFB$predicted[2] / dat_geoFB$predicted[1]
dat_geoFB$conf.low[2] / dat_geoFB$conf.low[1]
dat_geoFB$conf.high[2] / dat_geoFB$conf.high[1]

#Twitter Geographic references
dat_geoTW$predicted[2] / dat_geoTW$predicted[1]
dat_geoTW$conf.low[2] / dat_geoTW$conf.low[1]
dat_geoTW$conf.high[2] / dat_geoTW$conf.high[1]

#Facebook Regionalized wording
dat_regFB$predicted[2] / dat_regFB$predicted[1]
dat_regFB$conf.low[2] / dat_regFB$conf.low[1]
dat_regFB$conf.high[2] / dat_regFB$conf.high[1]

#Twitter Regionalized wording
dat_regTW$predicted[2] / dat_regTW$predicted[1]
dat_regTW$conf.low[2] / dat_regTW$conf.low[1]
dat_regTW$conf.high[2] / dat_regTW$conf.high[1]

